A method to accelerate computational efficiency by more than two orders of magnitude for Monte Carlo simulations of electron-solid interactions

A method has been developed to increase computational efficiency in Monte Carlo simulations of electron transport and interactions in matter. The method serves as the computational engine for the open-source code AMCSET (Aggie Monte Carlo Simulations of Electron and Ion Transport). The key is to combine n consecutive neighboring free flying distances into groups. Within each group, both flying distance and Mott scattering angles are obtained using Monte Carlo sampling under an equal energy approximation. This reduces the number of integrations of the tabulated differential Mott scattering cross-section in scattering angle selection, i.e., from 1000 to 1 if n = 1000. The method increases efficiency by more than 100 times. At the same time, the calculation still guarantees accuracy in calculating electron trajectory, excitation/ionization energy deposition, elastic scattering energy deposition, and displacement creation. For demonstration, 10 MeV electron bombardments of pure Fe with n up to 1000 are used as examples. The method, due to the availability of tabulated scattering cross-sections, is applicable for targets of the entire elemental table up to Z = 118, and for electron energies up to 900 MeV.

and energy after each collision is updated to reflect energy losses, including both continuous energy loss via inelastic scattering and discrete energy loss via Mott scattering.
Our method is applicable for electron irradiation but cannot be applied to other types of particle bombardments such as ion bombardment.This is due to the uniqueness of Mott scattering, which determines the trajectory of electrons, but whose energy loss is almost negligible compared to inelastic energy loss.This feature is important, allowing neighboring Mott scatterings to be evaluated under an approximately constant energy.If Mott scattering were a major energy loss mechanism, Mott scattering within the same group would become sensitive to its preceding scattering events, necessitating that scatterings be evaluated sequentially, one after another, as in the case of ion bombardment.Additionally, Mott scattering is not highly sensitive to energy changes, which allows for a constant energy approximation within a reasonably selected grouping size.Additional justification and details of the proposed method will be further explained in Sect."Modeling procedure".
The simulations in the present study consider ionization/excitation, bremsstrahlung, and Mott scattering.Ionization/excitation and bremsstrahlung are often treated as continuous energy loss, while Mott scattering is typically treated as discrete energy loss.The statistical fluctuations of the energy loss from excitation and ionization are ignored in the present study.To be precise, all interaction types can be treated as discrete energy losses and modeled as individual scattering events.In practice, free flying distances are introduced to simplify calculations by avoiding the need to compute each individual interaction event.A free flying distance represents the separation distance between two consecutive Mott scatterings.Within each free flying distance, energy loss due to excitation/ionization and X-ray production is calculated as the product of the flying distance and the corresponding stopping power.At the end of each free flight distance, a Monte Carlo approach is used to sample the deflection angle of Mott scattering.Below, the simulation details of each interaction mechanism are separately discussed.

Energy loss due to inelastic electron-electron scattering
For electrons with kinetic energies below 10 keV, corresponding to velocities significantly lower than the speed of light, the energy loss is calculated by a modified Bethe equation given below 17-19 : where.
N : the target atomic density in the unit of atoms/cm 3 .N a : the Avogadro's number ( N a = 6.02214076 × 10 23 ).Z : The atomic number of target atoms.E k : The kinetic energy of electron (in the unit of keV).J*: The modified mean ionization potential (in the unit of keV).s: the flying distance.
The modified mean ionization potential, J*, is calculated by 18 : For high-speed electrons, a relativistic form of Bethe equation is used 20 : where m e c 2 is the rest mass energy of an electron in the unit of keV ( m e c 2 = 511keV ).β is the speed of electrons normalized by the speed of light, β = ν/c.

Energy loss due to Bremsstrahlung
Bremsstrahlung is electromagnetic radiation which is significant when electron energies are significantly high.
The corresponding stopping power of electrons, measured in units of keV/cm, contributed by the Bremsstrahlung radiation, is calculated by 21 : (1) www.nature.com/scientificreports/

Mott scattering
For describing electron-nucleus scattering, the common approach is to modify the Rutherford scattering cross section with a Mott correction.The Rutherford scattering cross section is calculated by treating the target nucleus as an isolated point charge, without considering the presence of atomic electrons.The more precise Mott cross sections consider (1) the effect of electron spins during scattering, in addition to the Coulombic force, (2) relativistic effects for high-energy electrons, and (3) the electron screening effect.The latter is due to the presence of atomic electrons, which creates a shielding effect, reducing the effective charge of the nucleus to less than its atomic number.
Considering the aforementioned effects, the differential scattering cross section for electron-nucleus scattering can be described as follows 22 : where dσ Ruth d� is the Rutherford scattering cross section.R Mott is the Mott ratio which is for the correction after considering the electron spin effect but not electron-screening effect.The term 1 − F e q 2 is for the correction considering the electron screening effect.The parameter q is a function of electron velocity and Mott scattering angle, and is to be further explained shortly.An approximation of the Rutherford scattering cross section is given by 23 : where θ is the scattering angle of electrons, and r e is the classic radius of electron (r e = 2.817938 × 10 -13 cm).The Mott treatment is based on a method from Wenzel 24 , originally used to deal with incident and scattered waves on a point-like nucleus.After considering the spin effect, the obtained screening-free Mott differential cross sections are formulated as two conditionally convergent infinite series in terms of Legendre expansions 22 .A simple analytic expression of R Mott is available but it is limited to low Z materials 25 .For wider energy regions and more general target atoms, a fitting formula can be used with tabulated parameters.One generally accepted format is expressed in a polynomial expression 26 : β = 0.7181287 corresponds to the average velocity of electrons at 0.7181287c within the energy range of 1 keV to 900 MeV.α j and b k,j are fitting parameters which were originally proposed by Lijian et al. for elements up to Z = 90 26 .Boschini et al. provided an updated table for elements up to Z = 118 22 .
The term [1-Fe(q)] 2 , which accounts for screening effects, varies depending on the proximity of the electron to a target nucleus and the specific electron configurations of the target atoms.Consequently, it is dependent on both energy and the atomic number of the target atoms.If the Dirac-Hartree-Fock-Slater screening function is used, F e q can be expressed by 27 : where h is the Plank constant, A i and α i are parameters to describe the screening function, and q is the momentum transfer, calculated by 27 :

Determination of scattering angle under Mott scattering
Combining Eqs.(5-10), the total cross section for electron-nucleus scattering is calculated by the integration: Next, a random number R θ , which must be in the range of 0 to 1, is used to select the scattering angle at the end of each free flying distance, calculated by ( 5) where σ Mott (θ) represents the integration of the differential cross section up to the angle θ .In practice, the maximum allowable scattering angle, π, is divided into m intervals.The scattering angle θ falls in the region between θ i and θ i+1 if R θ satisfies where �θ i is the interval width for θ i .Since dσ Mott (θ ) d� is very high at small angles, θ i should be unevenly distributed, favoring a higher density at small angles.One such choice is: where m is the number of angle intervals used in the cross-section calculation, i is an integer in the range of 1 ≤ i ≤ m − 1 .�θ i for 2 ≤ i ≤ m − 2 is calculated as the distance between θ i−1 +θ i 2 and θ i +θ i+1 2 .At two boundaries, �θ 1 is the distance between 0 to θ 1 +θ 2 2 , and �θ m−1 is the distance between θ m−2 +θ m−1 2 and π.Once θ is determined to be in the region between θ i and θ i+1 through Eq. ( 13), its exact value needs to be more precisely determined based on the distance from the upper integration boundary θ i+1 +θ i+2 2 .It is important to note that this boundary value is slightly larger than θ i+1 due to the selection of �θ i+1 .The calculation assumes σ Mott (θ) and θ around θ i+1 follow the same linear proportionality.For 2 ≤ i ≤ m − 2, θ is calculated by: With Due to the high likelihood of R θ being less than σ Mott (θ1) σ T and the specific starting boundary selection for �θ 1 , the selection of θ between 0 and θ 1 is calculated by: with If θ falls into the last interval region, where σ Mott (θm−2) σ T < R θ ≤ 1 , then, due to the specific ending boundary selection, θ is calculated by: with The selection of θ for each Mott scattering represents the most time-consuming step in the entire Monte Carlo simulation.High accuracy necessitates calculating Eq. ( 12) with Δθ at an extremely small value (Δθ ≪ 1°).This requirement is due to the sharp changes of dσ Mott (θ ′) d� at very small angles.In Fig. 2a, dσ Mott (θ) d� of electron scattering by Fe at energies of 10 MeV, 1 MeV, and 100 keV is plotted.For 10 MeV electrons, when angles change from 0.1° degree to 1°, the differential cross sections are reduced by four orders of magnitude.This sharp reduction at extremely small angles results in a quick rise of σ Mott (θ) at small angles, as shown in Fig. 2b for θ < 1 • .For 10 MeV electrons, σ Mott (θ) quickly increases and saturates around a value of ~ 1.8 × 10 -18 cm 2 at θ > 0.2°.Normalized by the saturated value at large θ , σ Mott (θ) σ Mott (θ≫1) reaches about 0.5 at θ = 0.035 • and about 0.6 at θ = 0.042 • .This means that, to differentiate random numbers 0.5 and 0.6, the selection procedure requires θ to be differentiable at a resolution of 0.007 • .Equation (15) allows the extraction of θ even when it is less than the first interval width �θ 1 .However, this equation is based on the assumption of linear proportionality within �θ 1 .Therefore, �θ 1 cannot be wider than the critical θ value where the σ Mott (θ) value begins to deviate from this linear proportionality.The arrow in Fig. 2b marks the position of linear proportionality limit, which is 0.07 • .Based on Eq. ( 14) and the formula �θ 1 = 1 2 (θ 1 + θ 2 ), it can be determined that interval number n cannot be less than 1000.If an evenly divided interval is used, as opposed to the method in Eq. ( 14), the interval number increases to 2600.
This results in an incredibly large number of angle intervals for Monte Carlo simulations, significantly increasing the computation time.For lower energies, such as 1 MeV and 100 keV, the changes in dσ Mott (θ ′) d� at small θ are less pronounced.Additionally, the saturation values of σ Mott (θ) increases, and saturation begins at larger θ values.These factors relax the requirement of small �θ compared to the 10 MeV case.However, the number of angle intervals required remains significantly large, thereby slowing down the computation.For the 100 keV case, the σ Mott (θ) profile shows a shape different from both the 1 MeV and 10 MeV cases.Compared with 1 MeV, its proportional limit shifts to θ ≈ 0.012 • .Still, it requires a large n.

Determination of free flying distances under Mott scattering
The free flying distance of electrons must conform to Poisson distributions.Accordingly, the selection of the free flying distance, denoted as L, in Monte Carlo sampling for each individual Mott scattering event is calculated by: where σ T is the total Mott scattering cross section, as obtained from Eq. ( 11).R L is a random number created to determine free flying distances.Note that R L and R θ are different random numbers used for different purposes.
With N = 8.482 × 10 22 atoms/cm 3 for pure Fe, L is calculated for 10 MeV electron bombardment using Eq.(18). Figure 3a plots the distributions of L selections along the trajectory of one 10 MeV electron.Figures 3b-e plot the statistical distributions of L for electron energy ranges of 10-7 MeV, 7-4 MeV, 4-1 MeV, and less than 1 MeV, respectively.The distributions of L values become narrower and the mean free flying distance gradually reduces to short distances as the local electron energies decrease.This trend is attributed to the energy dependence of Mott cross section.The total cross section, σ T , increases with decreasing electron energy.The probability distribution profiles of Figs.3b-e follow the expected Poisson distribution functions.
Given that the total flying distance for a 10 MeV electron is approximately 0.8 cm (to be discussed shortly), the total number of L selections required for a single electron bombardment simulation is significant.The total number of free distance selections in Fig. 3a, for a single electron bombardment is approximately 90,000.For any meaningful Monte Carlo simulations, more than 1000 electron bombardments are likely required for reliable statistics.Additionally, considering that each free flying distance involves the integration of Mott scattering cross section with the number of angle intervals at 1000 or higher, the simulation becomes very time-consuming.This necessitates an innovative approach, which will be further explained in the next section.

Energy loss due to Mott scattering
To calculate the energy loss after each Mott scattering, the following general formula is used 28 : where M is the mass of the target atom.

Proposed new method to increase compactional efficiency
The key approach to improving computational efficiency is to reduce the number of times Mott scattering cross sections are integrated.As will be explained, this can be achieved by combining neighboring free flying distances into groups, assuming that each individual Mott scattering within the group can be simulated using the same electron energy.If so, only one integration of Mott scattering is needed for one group to determine all scattering angles.Below, we first provide justification for why such an approach is valid and appropriate.Then, the methodology of grouping free flying distances and calculating scattering angles within the group is introduced.Lastly, a comparison is made between traditional approaches and the proposed new methods to show the efficiency and accuracy of the new method.The equal-energy combination of multiple free flying distances is based on the fact that Mott scattering is most significant in changing the trajectories of electrons but much less significant in causing energy loss.www.nature.com/scientificreports/compares the stopping powers contributed by elastic scattering (Mott scattering), inelastic scattering (excitation and ionization), and bremsstrahlung, using electron bombardments in Fe as an example.The energy loss due to excitation and ionization is the most significant.The energy loss due to Mott scattering is systematically lower than that due to excitation/ionization by at least three orders of magnitude for the entire energy range.Bremsstrahlung becomes more significant than excitation and ionization at electron energies exceeding ~ 60 MeV.The insignificance of Mott scattering in energy transfer suggests that the local electron energy along their flying trajectory is not affected by the selection of Mott scattering angles.In other words, Mott scattering can be repeatedly evaluated assuming the same electron energy over a certain flying distance.As illustrated by Fig. 3b-e, the average free flying distance is on the order of 1 × 10 -6 cm.Hence, the equal energy assumption can be applied to, for instance, 1000 combined free flying distances.
One consequence of the almost negligible energy loss in Mott scattering (in comparison with the energy loss due to excitation and ionization, as shown in Fig. 4) is that the accumulated flying distances, integrated over the entire trajectory, are roughly the same for all electrons (as long as they are stopped inside the substrate).Figure 5a shows the trajectories of two 10 MeV electrons.Both exhibit large trajectory changes.Figure 5b plots the kinetic changes of both electrons as a function of flying distances.Black symbols and red symbols are used to differentiate the two electrons.Both electrons travel roughly the same distance of about 0.78 cm, although both have very large deflection histories as shown by the values in Fig. 5c. is the angle with respect to the z-axis.The same flying distances for both electrons are a consequence of the fact that kinetic energy loss is insensitive to Mott scattering.Even though individual scattering events show large differences in scattering direction changes and energy transfer, overall, Mott scattering does not cause significant differences in the total flying distances.
Another significant aspect of the results from Fig. 5b is the consistency in the slope, which appears to be similar throughout the entire flying journey.It's important to note that this slope represents the stopping power ( dE ds ).A consistent slope indicates that the stopping power for neighboring flight distance groups is roughly the same.This provides justification for why each Mott scattering event within a group can be calculated using the same energy level-specifically, the midpoint energy within the group.Although the actual energy within the group changes due to accumulated energy loss, such gradual changes do not significantly affect the stopping power.Hence, it is acceptable to use a singular stopping power value (a consequence of the equal energy approximation) to calculate the stopping loss across the entire group.
Our proposed method includes the following steps: (1) Define an integer number n as the size of the group, which is the number of flying distances included in each group.For the present study, n is up to 1000.(2) Produce n random numbers for each group.Each random number is in the range (0,1).
(3) Sort the random numbers in their original sequence from smallest to largest: R θ (i) < R θ (i + 1) for i from 1 to n − 1.  × n i=1 L i .(10) Subtract the kinetic energy of the electron by both the energy loss from Mott scattering and the energy loss from continuous non-Mott scattering.The final energy, E end(current) , is the new starting kinetic energy for the next group of flying distances: E end(current) = E starting(current) − �E Mott − �E non−Mott .(11) Calculate the ratio of the midpoint energy to the starting energy of the current group, using α = E starting(current) + E end(current) /2 /E starting(current) .(12) Use the α value to estimate the midpoint energy of the next group via E middle (next) = E end (current) × α .(13) Use the scattering angles from each Mott scattering within the current group to calculate the total accumulated deflection for the current grouped flying segments.The final scattering direction is the starting direction for the next group of flying distances.The calculation details will be discussed shortly.(14) Repeat the above cycles until the electrons either leave the substrate as backscattered electrons or stop inside the substrate.Stopping is determined based on whether the kinetic energy after the last collision is below a threshold energy for stopping, which is set to be 0.02 keV.
Figure 6 further explains step 4. For simplicity, only five random numbers are shown.These five randomly produced values of R within the range (0,1) are sorted to satisfy the relationship that R 1 < R 2 < R 3 < R 4 < R 5 .During the integration of dσ Mott (θ ) d� to obtain the solid curve in Fig. 6, corresponding θ i value for each R i are obtained.Since R i is in ascending, the integration process can identify all θ values at once.This is key to avoiding individual integration for each R. Afterwards, re-ordering of R within a group will be performed for the subsequent scattering simulations.
The schematics in Fig. 7 illustrate the definitions of angles.For step 13, after each Mott scattering within the group, the new directions, with respect to the z-axis of the Cartesian coordinate, and Φ with respect to the angle away from the x-z plane of the coordinate, are obtained using the following conversion: sin � cos � = cos� 0 cos� 0 cosϕ i sinθ i + cos� 0 sin� 0 cosθ i + sin� 0 sinθ i sinϕ i where 0 and 0 are the angles before scattering, and Φ and Θ are the new angles after scattering.θ i is obtained from step 4. ϕ i is a value randomly assigned in the range of 0 to 2π, based on the assumption that the substrate is amorphous.The above angle conversions will be repeated n times for each group.All energy losses are energy-dependent.The equal energy approximation within the same group requires an appropriate selection of electron energy.Using the energy midpoint of the group for scattering calculations is a straightforward choice.In the present study, an explicit method for such energy approximation is used, as described in steps 11 and 12.The energy midpoint of the next group is obtained using the final and starting energies of the current group.This method assumes that the ratio of the energy middle point to the starting sin � sin � = −sin� 0 cos� 0 cosϕ i sinθ i − sin� 0 sin� 0 cosθ i + cos� 0 sinθ i sinϕ i (20) cos� = cos� 0 cosθ i − sin� 0 sinθ i cosϕ i  www.nature.com/scientificreports/energy of the current group is the same for the next group.The simulation of the next group yields energy information that allows a new ratio to be calculated for the subsequent group.

Displacement production
The simulation easily integrates with calculations of various by-products of electron bombardment, such as displacement production.For that purpose, once the transferred energy to target atoms exceeds the threshold displacement energy, information about the primary knock-on, such as displacing direction and kinetic energy gained, is transferred to a subroutine.This subroutine is SRIM-like and capable of full damage cascade simulations.For the completion of the discussion, below the major steps involved in the full damage cascades are listed, and the results are shown in Fig. 10a-4, b-4, and c-4.
The ion-solid interaction subroutine has the following sequential steps for displacement calculations 29 : (1) A random number is produced to find the impact parameter using the following equation where N is the atomic density of the target atoms, R is a random number (2) Once the collision parameter b is obtained, it is used as input for the following equation to calculate the minimum separation distance between the knock-on atom and target atoms in the center of mass coordinate.
where the potential V is the ZBL potential 30 , which is widely used in ion-solid interaction community.It was obtained through fitting various ion-target systems, considering quantum mechanics effects.The ZBL potential is a binary potential that considers the electron screening effect between the ion and the target nucleus.This screening effect acts as a modification of the Coulombic interaction between ion and target nucleus.The ZBL potential is expressed as 30 : where Z 1 and Z 2 are the atomic numbers of the knock-on and the target nuclei, respectively, and r is the distance between the knock-on and target nuclei.The screening parameter a is described as 31 : where a 0 ( a 0 =0.529Angstrom) is the Bohr atomic radius.The screen function is calculated as: (3) The curvature of the radius at the minimum separation distance in the center of mass coordinates is calculated using the following equation: (4) The scattering angle of the knock-on in the center of mass coordinates is calculated using Biersack's "magic formula" 29 : where all parameters on the right are obtained from fitting and are calculable once ion-target systems and ion energies are known.
(5) The scattering angle is used to calculate the energy transfer to the target atoms, using the following equation 29 : where m 1 and m 2 are masses of the knock-on and target atoms, respectively.Transferred energy T will be subtracted from the Fe knock-on, accounting for nuclear stopping power loss.T is also used to judge whether additional target atoms will be displaced, depending on whether T is greater than the threshold displacement energy.
(6) The calculation is repeated for another flying distance, which is equal to the nearest neighboring distance.There is no need to use the free flying distance approach to save time, since the energy transfer from electrons to target atoms is very small.A 'monolayer' approach is sufficient and appropriate.The approach uses the average atomic distance of Fe as the flying distance between two scattering simulations.
Vol:.( 1234567890) www.nature.com/scientificreports/(7) Before the next ion-target scattering simulation, the electronic energy loss is calculated.This energy, which contributes to electron excitation in Fe, is very small for low-energy ion bombardment and is treated as continuous energy loss.More details about electronic energy loss calculation can be found in reference 29 .

Results and discussions
The proposed grouping of flying distances, even with n = 1000, does not cause noticeable errors in the energy loss of electrons along their flying paths.Figure 8a compares the electron kinetic energies for three individual bombardment events with n = 1, 300, and 1000.All three curves overlap well.Three boxes refer to the near-surface region, the middle, and the end of penetrations.For the near-surface region, as shown in Fig. 8b, all three curves exhibit the same amount of energy loss.The solid line, hollow circle, and solid square refer to n = 1, 300, and 1000, respectively.In the middle range regions, as shown in Fig. 8c, all three curves still follow each other closely.At the end of the projected ranges, as shown in Fig. 8d, they are still indistinguishable.This is a consequence of Mott scattering not being a significant energy loss mechanism.Hence, individual electron trajectories, even with large variations in Mott scattering details, do not bring noticeable differences to the total flying distance.If Mott scattering details are not able to differentiate one trajectory from another, it justifies our grouping method by ignoring the differences between neighboring flying distances and by using the approximated same energy within the group for scattering angle selection.Note that all three electrons stop their penetrations at the same distance when the energies are below 20 eV.
The grouping method does not cause noticeable errors in the scattering angle statistics distributions either.Figure 9a plots the θ distributions for 1000 electron bombardments as a function of their flying distance, for the case of n = 1 (no grouping).The color indicates the density of data points, increasing from blue to red.θ values gradually increase from 0, as the initial bombarding conditions, to larger numbers with increasing penetration.At the very end of the penetration, θ becomes widely distributed over the whole 0 to π region.Figure 9b is a comparison with Fig. 9a, but for the case of n = 1000.Since flying distances are grouped, scattering analysis occurs after every 1000, which leads to a discrete pattern, instead of random distributions.At the very end, a similar spreading over 0 to π was observed.Since the free flying distance for each Mott scattering becomes smaller at low energies, the distance between neighboring groups becomes smaller, and the data density increases.The A comprehensive comparison is made in Fig. 10 to show the equivalence and accuracy of the L grouping method for the calculations of electron distributions, energy deposition, and displacement production.Figure 10a-1 is the plot of the electrons as a function of their radius and penetration depth in polar coordinates for the case of n = 1 (no grouping).We selected polar coordinates since the yield is symmetrical.In a comparison, similar 3D distribution patterns are obtained for n = 300 (as shown in Fig. 10b-1) and for n = 1000 (as shown in Fig. 10c-1).Both the pattern shape and peak location agree for the three cases.The patterns are characteristic of particles which experience backscattering.Figure 10a-2, 10b-2, and 10c-2 compare the deposition of ionization/ excitation energy for n = 1, 300, and 1000, respectively.All three agree with each other, featuring the highest energy depositions at the very beginning of the bombardment.Figure 10a-3, 10b-3, and 10c-3 show energy deposition for elastic scattering (Mott scattering).All patterns agree with each other and are very similar to the pattern of inelastic scattering.As discussed already from Fig. 4, elastic scattering energy loss is more likely a systematic shift of inelastic scattering energy loss as a function of energy.Therefore, it is no surprise that both share similar patterns, but the amplitude of elastic scattering energy loss is systematically lower than that of inelastic scattering energy loss.Figure 10a-4  The displacement production is a byproduct of Mott scattering, and occurs only when the transferred energy to target atoms exceeds the displacement threshold energy (40 eV for Fe).Hence, the displacement pattern is similar to that of elastic energy loss pattern but has a cutoff when the energy loss density is below a threshold, which leads to 'shrunken' patterns in comparison with elastic energy deposition patterns.
Figure 11a-d compares the one-dimensional distributions of electron distributions, excitation/ionization energy deposition, elastic scattering energy deposition, and target atom displacement production electrons as a function of depth for 10 MeV electron bombardments in Fe.The profiles were obtained using the 3D profiles in Fig. 10.All these plots agree well for three cases: n = 1, 300, and 1000.
Figure 12 summarizes the computation time required for the simulations, taking flying distance group sizes as n = 1, 10, 100, 300, and 1000.All times are normalized against that taken for n = 1000.There is a significant improvement in efficiency when n changes from 1 to 100.Efficiency appears to saturate at n ≥ 300.Comparing n = 1 with n = 1000, it is evident that the proposed method increases efficiency by approximately 400 times.All calculations were obtained using an angle interval number of 1000.
In the present study, the Mott differential scattering cross-section is calculated using R Mott obtained from fitting 26 .The fitting formula with tabulated parameters is very helpful in reducing computation time and broadening the applicability for a wider energy region and more general target atoms.The error is less than 1.0% for atomic numbers Z = 1 to 44 for energies from 1 keV to 900 MeV.
To assess the differences between the present study and the exact solutions from other groups, a comparison is made in Fig. 13.The solid line shows the results combining fitting-obtained R Mott and the Dirac-Hartree-Fock-Slater screening function 26,27 , and the dashed line shows the exact solution calculated from the Dirac-Hartree-Fock potential for 10 keV electron irradiation in nickel 31 .The solid line is slightly higher than the dashed line at angles less than 70 degrees.For angles greater than 70 degrees, the dashed line becomes slightly higher than the solid line.The deviation of the two curves over the entire range is less than 40%.Figure 13 shows that for the most important scattering angle range, less than 5 degrees, the two curves are almost identical.The color refers to the log scale of the corresponding values.For electron distributions, the value is the number of electrons per cm 3 per incident electron.For inelastic and elastic energy deposition, the values are in the unit of keV per cm 3 per incident electron.For the Fe displacement plot, the unit is the number of Fe displacements per cm 3 per incident electron.All plots are obtained using 3000 electron bombardments.The number of angle intervals is 1000.

Conclusion
We proposed a new method to speed up Monte Carlo simulation of electron-solid interaction which considers inelastic scattering, elastic scattering, and bremsstrahlung, for energy ranges up to 900 MeV.The key approach is to maximize the accuracy of integration of differential Mott scattering using small angle intervals and to minimize the number of integrations.The latter is a key limiting factor in computational speed.The method combines free flying distances into groups, with each group containing up to 1000 free flying distances.The scattering angle selections for each free flying distance within the group are achieved using only one integration process.The grouping method assumes that electron energy within each group is approximately constant.This assumption is supported by the fact that Mott scattering contributes an energy loss which is a few orders of magnitude lower than that from excitation and ionization.The grouping method does not affect the total flying distance nor the statistical distribution of scattering angles.For the extreme case of grouping 1000 neighboring free flying distance for each group, the computation speed increases by more than two orders of magnitude, while the obtained 3D distributions of electrons, energy depositions, and displacement creations still agree well with the simulation without grouping.

Figure 1 .
Figure 1.Schematics of scattering calculation used (a) in the traditional approaches and (b) in the present study.

Figure 2 .
Figure 2. (a) Differential cross sections as a function of scattering angle θ in Fe and (b) integrated cross sections from angle 0° to θ, for 10 MeV (dark solid line), 1 MeV (blue dashed line), and 100 keV (red dotted line) electrons.

Figure 3 .
Figure 3. (a) The selection of free flying distance as a function of local electron energy for the simulation of one single electron bombarding pure Fe, at an incident energy of 10 MeV.(b-e) the statistics analysis of L selection for the electron energies in the range of 10-7 MeV, 7-4 MeV, 4-1 MeV, and < 1 MeV.The number of angle intervals is 1000.

Figure 4 .( 4 )
Figure 4.The stopping power (energy loss per unit flying distance) contributed by inelastic scattering (excitation and ionization), bremsstrahlung, and elastic scattering (Mott scattering) in Fe.

Figure 5 .
Figure 5. (a) 3D trajectories of two electrons in Fe, both at an incident energy of 10 MeV; (b) local kinetic energies as a function of flying distances for both electrons; (c) local θ values with respect to the z-axis as a function of flying distance for both electrons.Each line in (b) contains about 90,000 points (free flying distances).The number of angle intervals is 1000.

Figure 6 .
Figure 6.Obtaining θ i values for each R i for a group of five random numbers in ordered sequence from low to high.

Figure 8 .
Figure 8.(a) Electron kinetic energy as a function of flying distances for three group numbers, n = 1, n = 300, and n = 1000.(b) Enlarged region near the surface, (c) enlarged region in the middle of the penetration, and (d) enlarged region at the end of their ranges.The number of angle intervals is 1000.

Figure 9 .
Figure 9. (a) Density plot of scattering angles for 1000 electron bombardments in Fe as a function of their flying distances, for the case of n = 1 (without L grouping), (b) density plot of scattering angles for 1000 electron bombardments for the case of n = 1000, (c) statistics analysis of scattering angles for both n = 1 and n = 1000 over eight distance regions of 0-0.1 cm, 0.1-0.2cm, 0.2-0.3cm, 0.3-0.4cm, 0.4-0.5 cm, 0.5-0.6 cm, 0.6-0.7 cm, and > 0.7 cm, respectively.The number of angle intervals is 1000.

Figure 10 .
Figure 10.(a-1, a-2, a-3, a-4) 3D plots of electron distributions, excitation/ionization energy deposition, elastic scattering energy deposition, and target atom displacement production for 10 MeV electron bombardment in Fe, for n = 1 (no grouping of flying distances), (b-1, b-2, b-3, b-4) plots of the same but for n = 300, (c-1, c-2, c-3, c-4) plots of the same but for n = 1000.The color, from blue to red, indicates the value change from low to high.The color refers to the log scale of the corresponding values.For electron distributions, the value is the number of electrons per cm 3 per incident electron.For inelastic and elastic energy deposition, the values are in the unit of keV per cm 3 per incident electron.For the Fe displacement plot, the unit is the number of Fe displacements per cm 3 per incident electron.All plots are obtained using 3000 electron bombardments.The number of angle intervals is 1000.

Figure 11 .Figure 12 .
Figure 11.One-dimensional depth profiles of (a) electron distributions, (b) excitation/ionization energy deposition, (c) elastic scattering energy deposition, and (d) target atom displacement production for 10 MeV electron bombardment in Fe, for n = 1, 300, and 1000.All results are normalized to one electron.

Figure 13 .
Figure 13.Comparison of the Mott scattering differential cross sections obtained from the study and Reference 31 .